A GENERAL APPROACH TO ENHANCE SLOPE LIMITERS ON 

NON-UNIFORM GRIDS 
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Abstract. A general approach to study and enhance the slope limiter functions on non-uniform 
grids is presented. Slope limiters are preferred in high-resolutions schemes in general and MUSCL in 
particular to solve hyperbolic conservation laws. However, most ID limiters are developed assuming 
uniform meshes in space, which are shown to be inadequate on non-uniform grids. Especially, second- 
order convergence is shown to be lost when the conventional limiters are applied on irregular grids 
in the case of smooth solutions. A methodology based on the classical reconstruct-evolve-project 
approach and Harten's stability theory is presented to study the slope limiters on ID non-uniform 
computational grids. Sufficient conditions for the limiters to lead to formal second-order spatial 
accuracy, total- variational-diminishing stability and symmetry-preserving property are derived. The 
analysis and results extend naturally to cell-centered finite volume methods in multiple dimensions. 
Several most widely used conventional limiters are enhanced to satisfy these conditions, and their 
performances are illustrated by various ID and 2D numerical examples. 
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1. Introduction. Limiting strategies are often employed in high-resolution finite 
volume methods (FVM) to solve scientific and engineering problems that are governed 
by time-dependent hyperbolic conservation laws (HCL). These problems are usually 
characterized by discontinuous solutions at finite time even though the initial condi- 
tions are smooth. In practice, two classes of limiting techniques are often employed: 
flux limiting, exampled by Flux-Corrected Transport (FCT) methods [5J [371 IS] ; and 
slope limiting, like the Monotone Upstream-centered Schemes for Conservation Laws 
(MUSCL) schemes [23l [131 HI [10] . Both techniques try to keep the numerical scheme 
high-order accurate (often second-order) everywhere in the computational domain ex- 
cept near discontinuities, where order of accuracy is lowered to first-order in space to 
achieve certain stability property, like total variational diminishing (TVD); and they 
are equivalent to each other in certain situations. This paper focuses on slope limiters, 
but the latter are also reformulated in the analysis part of the paper as equivalent 
flux limiters to provide better insight of their functionalities. 

Despite their wide acceptance in practical computation, the most widely used 
slope limiters are developed in the context of ID uniform meshes. These limiters 
are applied to solve problems in multiple dimensions by: (1) dimension-splitting on 
structured meshes, or (2) solving a local directional ID problem on unstructured 
meshes. On the other hand, the limitations of these conventional limiters are well 
documented in literature [24l [^ ■ V. Venkatakrishnan [24] noticed that non-smooth 
limiter functions led to non-converging solutions when considering steady problems. 
Thus smoothness of the limiter functions is desired in practical computations. M. 
Berger et. al. [2] reformulated the slope limiters in a symmetric form other than 
the classical one, and studied their effects on non-uniform grids. In this report, 
conventional limiters are shown to lead to less than second-order accuracy on non- 
uniform grids; and several modifled limiter functions, including two modified van Leer 
limiters, are provided. However, the proposed generalized van Leer limiters therein 
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lose the C°° property that is possessed by the conventional one; in addition, no 
generalization to the popular van Albada limiter is provided. 

In this paper, a different approach is taken to analyze the functionality of limiter 
functions when non-uniform grids are employed in computation. This approach mod- 
ifies the original reconstruct-evolve-project (REP) procedure [7] to take into account 
of the irregularity of the computational grids; and it requires no reformulation of the 
variables (cf. [1]). The analysis focuses applying MUSCL scheme to ID problems on 
non-uniform grids, in which case sufhcient conditions are derived for designing the 
limiter functions such that TVD stability and formal second-order accuracy in space 
are achieved. In addition, symmetry-preserving property on non-uniform meshes are 
formally defined, and a sufficient condition for the limiters to preserve symmetric 
solutions is provided. The analysis in ID extends naturally to multiple dimensions 
for cell-centered FVMs [TJ [5] . The reason that it is not suitable for vertex-centered 
FVMs, for example, the dual-cell approach [T], is that these methods usually compute 
fluxes at edge centers, which is equivalent to a locally uniform ID grid. 

Several most widely used slope limiters are enhanced to satisfy the derived con- 
ditions. The enhancements include novel modifications to the van Leer limiter and 
the van Albada limiter, which preserve the C°° property of their conventional coun- 
terparts. The performance of these enhanced limiters are illustrated by various ID 
and 2D numerical examples. 

The remainder of the paper is organized as follows. In Section [2] FVM with 
MUSCL is briefly reviewed and a case study of solving ID Euler equations on non- 
uniform grids is presented to raise the issue of concern. Mathematical analysis of 
applying slope limiters with non-uniform ID meshes are given in Section [3j where 
sufficient conditions are derived for designing slope limiter functions such that the 
resulting schemes are TVD stable, second-order accurate in space and preserve sym- 
metric solutions. Several most widely used limiter functions are enhanced to satisfy 
these conditions, and they are given in Section |4] This section also shows by a case 
study of the van Leer limiter, which explains why present research is important. The 
performance of the enhanced limiters is evaluated by numerical examples in Section [5j 
Finally Section [6] concludes the paper. 

2. Case study. 

2.1. The MUSCL scheme. Consider the MUSCL scheme [23] that is applied 
to solve the general ID scalar conservation law 

ut + /(m)x = (2.1) 

In typical FVM, the ID computational domain 17 is divided into non-overlapping 
cells ili = [xi_i/2,Xi+i/2], where ••• < Xi_i/2 < 2:^+1/2 < Xi+3/2 < ■■■ are cell 
faces. The cell center of is defined hy Xi = ^ (2^1-1/2 + 2^^+1/2), and Axi := \ \ = 
Xi+1/2 ~ Xi_i/2 is the size of the cell. The mesh is supposed to be fixed in time. 

In the remaining of the paper, the following notations are used: u{x, t) denotes 
the exact solution of the PDF of a generic variable m; Ui{t) designates the exact 
cell-averaged data at time- instance t over the cell 17^; m" designates the numerical 
approximation of the value Ui{t"'); and finally Ui denotes the semi-discretized approx- 
imation of Ui. 

Integrating Fq. ( |2.1[ ) in space over fli one obtains 
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where Ui{t) is defined by 



1 



u{x, t)dx 



(2.3) 



Eq. (2.2) leads to the semi-discretized equations for the approximated solution u, 





duj 
~dt 



i+l/2 



F, 



i-l/2 



(2.4) 



Here is the numerical flux approximating f{u{xi^i/2,t)); and it is calculated 

using the semi-discrete solutions Ufc. The canonical Godunov scheme [6 defines the 
numerical flux -P^-f 1/2 as the solution of the Riemann problem defined by left data Ui 
and right data u^+i 



TT^Godunov TpBiem 

^i+1/2 - ^ 



{Ui,Ui+i) 



Whereas due to the high computational cost of the exact Riemann solutions, F^*^™ 
is often replaced by more economical approximate Riemann solvers, for example, the 
Roe solver [H] F^"". 

The spatial approximation i^j;-|-i/2 = F^°^{ui,Ui-^-l) leads to first-order error in 
space, and a second-order extension is to linearly reconstruct the state vector to the 
left and to the right of the cell face Xi+1/2 



F, 



+1/2 = F^°''{ui + ^crjAa;,,,u,,+i - ^di+iAxi+i) 



(2.5) 



where ai^i+i are the slopes for the linear reconstruction in ili^i+i. In order that 
Eq. (2.5) gives rise to a second-order accurate discretization in space, Oi should be 
consistent with the local gradients Ux at the cell centers Xi [13], denoted by D^Ui. 
One of the most natural choices for at on uniform meshes is given by 



(7,: = D^Ui -^(ui+i - Ui) 
Ax, 



(2.6) 



To avoid spurious oscillations near discontinuities caused by Eq. (2.6), the MUSCL 



approach multiplies D^Ui with a restricting factor </>, called slope limiter 



(2.7) 



Here and in the remains of the paper, the symbol <f> is reserved for slope limiters. 

Once the limited slopes Ui are computed, one use any preferred ODE solver, like 
forward Euler (for simplicity of analysis), to discretize Eq. (2.4) in time 



i+1/2 ^i-1/2 



At 







(2.; 



Integrating Eq. (2.2) over the time interval ] one has the exact integral 



form of ( 2.1 ) 



e,(t"+i)-M,(0 
At" 




f{u{Xi+i/2,t))dt 



I{u{x,^l,2.t))dt\ =0 



(2.9) 
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where Ai" = — t". It is clear that the numerical flux ^7^i/2 represents an 
approximation to the time-averaged flux at the cell face 



i+l/2 



1 



(2.10) 



2.2. Slope limiters. In practice, 0^ in Eq. (2.6) is a function of a local smooth- 



ness monitor 9i , of which a popular choice [5D1 12] is given by 

t>i = 

Ui+1 - Ui 



(2.11) 



this is the default smoothness monitors throughout the paper unless otherwise men- 
tioned. Conventional practice of MUSCL defines the limiters (/),; as 



(2.12) 



And the absence of the subscript of (p in the right hand side of Eq. (2.12) highlights 



the common strategy that a uniform limiter function is applied everywhere in the 
computational domain. A. Harten 9J showed on uniform meshes, a sufficient condition 



for the fully discretized system Eq. (2.81 to be TVD stable is 



< 



< 2, < 



< 2, V0 e 



(2.13) 



On the other hand, the effectiveness of limiters in retaining formal high-order 



accuracy when the solutions are smooth can be studied by rewriting Eq. (2.5) in the 
flux-corrected form 



F{ui + ^aiAxi,Ui+i 



/2 



-1/2 (^/^l 



/2 



-1/2 



(2.14) 



Here is any first-order numerical flux -Fi^x/2 ~ F{ui, u^+i), obtained from an exact 
or approximated Riemann solver; and F^ designates a high-order flux. The multiplier 
(p defined at the cell face 2:^+1/2 is thusly a switch between the high-order flux and 
the low-order flux, and hence called flux limiter. Typically, ipi+i/2 is a function of 
another local smoothness monitor defined at cell faces 9i+i/2 



(2.15) 



where the absence of the subscript of ip in the right hand side, again, highlights the 
practice that the uniform flux limiter is employed everywhere. 

In order to preserve high-order spatial accuracy, a sufficient condition is ipi+1/2 = 
1 for linear data (linearity preserving); in the case of uniform grids, this condition 
often reduces to 



^(1) = 1 



(2.16) 



It can be shown for linear problems that g iven the smoothness moni tor 6 j, one 
can define ^i+1/2 and F^ properly so that Eq. (2.141 is equivalent to Eq. (2.5|. This 
equivalence provides a convenient tool to study the formal order of accuracy a limiter 
could deliver, as presented in Section [3] 
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Several the most widely used limiter functions are listed as follows 
minmod: 



LTmmnod 



superbee: 
MC: 

van Leer: 
van Alhada: 



°'^(6') =niax(0,min(6',l)) 
^superbee ^Q-^ = max(0, niin(26', l),min(6l,2)) 

b^'^{e) = max(0, min(26', ^(1 + 61), 2)) 

+ 1^1 



ivan Leer 



(0) 



ivan Albada 



(0) 



1 



(2.17a) 
(2.17b) 

(2.17c) 
(2.17d) 

(2.17e) 



All the five limiter functions satisfy the conditions given by Eq. (2.13) and Eq. (2.16). 
Among the five, the minmod [TS], superbee [TS], monotonized central [MC] [22] are 
piecewise linear; and the van Leer limiter |21j and the van Albada limiter [19j are 
smooth functions of the smoothness monitor. 

2.3. Loss of accuracy. The case study is to solve the ID Euler equations with 
periodic boundary conditions 



Wt + f(w)a; =0, W 



P 

pu 
E 



f(w) = 



pu 
pu^ + p 
u{E + p) 



xe[-l,l] (2.18) 



w is the conservative fluid state vector and f is the physical flux function, p, u and p 
are density, velocity and pressure respectively; E is the total energy per unit volume 



E ^ pe+ -pu 



P 



(7-l)p 



7 



1.4 



Here e is the speciflc internal energy and is given by the ideal gas law. The initial 
condition is the following and also plotted in Figure |2.1[ 



(2.19) 



p{x, 0) = 1 + 0.5 sin(7rx) 
u{x,0) = 2 + 0.5sin(7rx) 
p{x, 0) = 1 + 0.5 sin(7ra;) 



Both uniform and non-uniform meshes are tested. The uniform mesh has N equal- 
size cells: —1 = Xi/2 < x^/2 < ■ ■■ < X(2N-i)/2 < ^{2N+i)/2 = 1- Thc non-imiform 
mesh is constructed by perturbing the uniform grid as follows 

1. Given N, let x'^^-^^^i * = 0, 1, • • • , iV be the cell faces of the uniform mesh. 

2. Set a fixed ratio r : < r < 0.5. 

3. For each inner grid point 2:^^1/2'* — 1; ' ' ' i-^ ^ 1: define Xi_i_i^2 — ^[+1/2 + 
rdi^i^2 where 6^^1/2 are independent random variables with uniform distri- 
bution on [-h, h], h = 2/N. 



The two cndpoints are fixed: Xi/2 — ^1/2 ^^"^ 2^(2^+1)72 



cii^^ ■^{2N + l)/-2 - ^{2N+l)/2- 

MUSCL with Roe flux and limiter functions in Eqs. (2.17) in space, and second- 
order TVD Runge-Kutta (TVD RK2) 8J in time arc employed in all tests. A CFL 
number 0.6 is used with the reference cell size 2/N for both uniform and non-uniform 
meshes. This relatively small CFL number is chosen since cell sizes in the non-uniform 
grids can be smaller than the reference size. 
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Fig. 2.1. Initial condition for the test problem FiG. 2.2. Convergence rates in pressure 

achieved using van Albada limiters on various 
meshes 

Table 2.1 

Convergence rates achieved on uniform meshes by applying conventional limiters 





minmod 


supcrbce 


MC 


van Leer 


van Albada 


p 


1.6671 


1.9578 


2.0036 


2.1405 


2.2246 


u 


1.7368 


1.6901 


2.0703 


2.2474 


2.2716 


p 


1.7235 


1.8475 


2.0067 


2.1725 


2.3348 



Given any limiter function and perturbation ratio r, solutions on five grids with 
sizes ranging from 100 to 1600 are calculated. The global Li errors in pressure cal- 
culated with conventional van Albada limiter are reported in logarithmic scale in 
Figure |2.2[ and the convergence rates for density, velocity and pressure are reported 



in Tables. [2^1 2.3 The details of calculating the convergence rates, especially on the 
non- uniform meshes, are offered in Appendix [A] 

Clearly, second-order convergence is lost when the non-uniform meshes are em- 
ployed; and the obtained convergence rates are degraded to the first-order. On the 
other hand, none of the simulations fails (up to T = 6.0, see Figure 2.3), indicating 
stability is still retained. These are exactly the issues this paper concerns. 

Remark: practical implementations sometimes involve a switch mechanism that 
turns on the slope limiter only if discontinuities are detected. Nevertheless, such 
a strategy falls into the scope of this paper for the following reason. Consider the 
switch depends on some smoothness monitor 9, which satisfies = 1 in the case of 
linear data. Suppose the switch is turned on if |0 — 1| > 6 for some small positive 
number S; then given any preferred slope limiter function this switching strategy 
is equivalent to applying the following limiter 



^switch {S) = X{\0-i\>5}<t>{S) +X{\0-i\<s} 

where x is the characteristic function. The issues raised in this section equally applies 
to the slope limiter function (j) switch- 

3. Mathematical analysis. 

3.1. Slope limiters and TVD stability. The classical REP approach is ap- 
plied to study the effectiveness of slope limiters on non-uniform ID grids to solve the 
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Table 2.2 

Convergence rates achieved on non-uniform meshes with r = 0.2 by applying conventional limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


p 


1.0993 


1.0562 


1.0918 


1.1426 


1.1401 


u 


1.1186 


1.0629 


1.0861 


1.1425 


1.1764 


p 


1.1051 


1.0407 


1.0652 


1.0915 


1.1124 



Table 2.3 

Convergence rates achieved on non-uniform meshes with r = 0.3 by applying conventional limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


p 


1.1015 


1.0895 


1.1190 


1.1445 


1.1469 


u 


1.1240 


1.0060 


1.1293 


1.1850 


1.1829 


p 


1.1084 


1.0215 


1.0812 


1.1215 


1.1331 



model problem given by the ID scalar advection equation 

Ut + CUx = 



(3.1) 



where c is a constant. Assuming non- uniform mesh with cell centers denoted by Xi as 
before, the REP procedure for updating the data from to is as follows 

1. The cell-averaged data u" correspond to from previous time step. 

2. (Reconstruct.) The data are linearly reconstructed on each cell to give the 



function u{x,t") (as illustrated in Figure 3.1 and Figure 3.2) 

u(a;,r) = u" + CTf (x - Xj), x e {xi^i/2, Xi+i/2) 



(3.2) 



3. (Evolve.) Eq. (3.1) is solved exactly over the time interval [t",t"+^] with 
initial condition given by Eq. (3.2 1; and the solution is denoted by u(a;,i"+^) 
(as illustrated in Figure 3.3). 

4. (Project.) Cell-averaged values at time t"+^ (Figure 3.4) are calculated as 



1 

Ax 



u{x,e+')dx 



(3.3) 



For general HCL, the "evolve" and "project" steps may be performed approximately. 
A classical result concerning the TVD stability is given by Harten: 
Theorem 3.1 (Harten). // the numerical scheme can be written, for one time 

step, in the following form: 



Then a sufficient condition for the scheme to be TVD is: 

Cr„i>0 yi,n 
D'l > Vi, n 
C" + £>" < 1 Vi, n 



(3.4) 



(3.5) 



here Cf, Z?" can be any numbers including those that are data dependent. 

In his original derivations, Harten assumed uniform grids; however this theorem can be 
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-O.S -0.6 -0.4 -0.2 0.2 0.4 0.6 0.3 1 



Fig. 2.3. Densities at T = 6.0 using minmod limiter (left) and van Albada limiter (right) on 
uniform and non-uniform meshes 



Fig. 3.1. REP-1: Cell-averaged data at t" 



Fig. 3.2. REP-2: Linear reconstruction 



applied to non-uniform meshes naturally, by observing that the definition of discrete 
total variance (Eq. (|3.6|) is independent of the particular cell sizes. 



(3.6) 



thus the original proof of Theorem 3.1 carries directly to any ID grids 



Next, the exact solutions to Eq. (3.1) obtained from Eqs. (3.2)-(3.3) are 



,n+l 



1 



= < - - - 7^^^ [(1 - \)Ax,a'^ - (1 - A,_i)Aa;,_iar_i] (3.7a) 



if c > 0, and 



,n+l 



- A.« - <+i) + -A, [(1 - AOAx.ar - (1 - A,+i)Aa;.+ifTr+i] (3.7b) 



if c < 0. Here Xi = |cAt"|/A2;i is the absolute value of the local Courant number. It 



is supposed here that < A; < 1. Assuming the smoothness monitor Eq. (2.11) and 



the approximated gradients Eq. (2.6 1, the limited slopes are 



1 



(3.8) 



Note that the limiter function (j) is equipped with subscript i to allow variance among 



cells. Plugging Eq. (3.8 1 into the exact solutions given by Eqs. (3.7), one can choose 

as 



the parameters C" and ZJf in Theorem 3.1 



Cr_i = A, + ^A.(l - A,)^ - ^A,(l - A._i)0._i(&,_i) 



D? = 



(3.9) 
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if c > 0, and: 



,)--A,(l-A,+i)- 



(3.10) 



if c < 0. Then the theorem leads to the conclusion that the REP approach is TVD 
stable given that 

0<Cr_i<l, ifoO 
and if c < 



which are satisfied given that 



Q<4>i{e)<2, 0< 



<2, Vi,( 



(3.11) 



and < A, < l,Vi. 

Indeed, in the case c > 0, assuming the above conditions, one has 

Cr_i > A, + iA,(l - AO • - iA,(l - A,_i) • 2 = A,A,_i > 



and 



Cti < A. + ^A,(l - A,) • 2 - iA,(l - A,_i) • = 2A, - A? < 1 



which indicates TVD stability by Theorem 3.1 The proof for < < 1 in the case 
c < is similar. 

To conclude, it is shown that a necessary condition for the limiter functions to 
lead to TVD stability on non-uniform grids is given by Eqs. (3.11), which are the 
same as Eqs. (2.13). It thusly explains why TVD stability is retained, even if the 
conventional slope limiters are applied with non-uniform meshes, as seen in Section [2| 

3.2. Second-order spatial accuracy. To investigate the formal spatial order 
of accuracy of the REP approach, the following fluxes are introduced 



■i+l/2 



c< + \c{l - A,)Aa:,crf c> 

c<+i - ^c(l - A,;+i)Aa;i+icrf_^i c < 



(3.12) 



so that the exact solutions given by Eqs. (3.7) is equivalent to those obtained from 
the finite volume formulation 
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Plugging the limited slopes (3.8 1 into Eq. (3.12) leads to 



-1/2 — 

when c > 0, and 



+ fc< + \b,+„2c{1 - A,)«+i - <) - c<) (3.13) 



i+l/2 



i+i H 5 cw,+i - -Bi+i/2c(l - Aj+i)(Uj+2 - u^+i) - C"»+l 

-Di+1/2 V ^ 

(3.14) 



when c < 0. Here -Bi+1/2 are coefficients to be determined. 



In orde r to express the fluxes given by Eqs. ( 3.13 )-( 3.14) in the flux-corrected 
form (2.14), the low-order flux and high-order flux need to be defined. In 
particular, F^ is chosen to be the first-order upwind flux 



cu- 



F. 



+1/2 



c> 
c < 



(3.15) 



and accordingly the high-order fluxes are 

:< + ii?.+i/2C (1 - «+i - <) 



F- 



+ 1/2 



C<+i - iS.+ l/2C (1 + ^) «+2 - <+l) 



c> 



c < 



(3.16) 



To this end, -Bi+1/2 should be chosen such that the fluxes (3.161 are second-order 
accurate in space. In the view of Eq. (2.10), and by applying Taylor series expansion 
in time to evaluate the time-averaged flux in the right hand side of Eq. (2.101 



1 

At" 



Cu{Xj+i/2,t)dt ■ 



1 

At" 



j„4 



c [u* + {t~ t"X] dt + oiiAf"y) 



CM* + -cAt"< -I- 0((At")2) 



(3.17) 



Here the superscript * is used to indicate that the values are evaluated at t"). 
By feeding the exact data at t" to the fluxes in Eqs. ( |3.16 1, and noticing that u{xi, t) = 
u{xi,t) + 0{Axf) and u{xi^i,t) = u{xi+i,t) + 0{Axf^i), one obtains 



B 



i+l/2 



0{Ax',+AxD 



if c > and 



T~l J. J. .It. "^1 

i-l-1/2 — + 2*^'^^ 



,-Bj+l/2 1 + 



cAt" 
Ax, 



cAt" 
Axj+i 



(Axj -I- Axi+i) - Axi 



(3.18a) 



(Axj+i + Ax.,+2) + Ax^+i 



OiAx^+^ + Ax-^^2) 



(3.18b) 



if c < 0. Comparing Eqs. ( 3.18 ) and Eq. ( 3.17 ), and applying -I- cu* = 0, the 
desired value i?i+i/2 is 

2Axi 



Bi+1/2 



Axi+Axi-f-i 

2Axi^i 
Axi^i+Axi^2 




c < 



(3.19) 
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which leads to the following flux limiters comparing Eq. (2.14) and Eqs. (3.13H(3.14| 

V^+l/2 - <^ (3.20) 

c < 



2A 



On the uniform meshes, the equation above confirms that the slope limiters and flux 
limiters are equivalent to each other. This equivalence is, however, not true on non- 



uniform meshes. Given linear data, one has from Eq. (2.111 

- Ui^i Ax,i_i + Axi 

Ui — ~ — 

Ui+i - Ui Axi + Axi+i 
Hence a sufficient condition for second-order space accuracy is <^i+i/2 = 1 or 

' Axi^i + Axi\ 2Axi 



Axi + Ax 



Axi + Axi+i ' 



(3.21) 



Thus in appearance of non-uniform grids, limiter function (f)i{-) must be defined locally, 
such that local grid geometry is taken into account. 



Finally, the TVD stability condition given by Eqs. (3.11 1 and the order condition 



given by Eq. (3.21 ) are always compatible with each other as shown below 

2Axi 



< 



^ 2Axi + 2Ax^+i _ 2 



Axi + Axi+i Axi + Axi+i 



and 



P Axj-i + Aa;^ ^ 



2Aa 



Axi + Ax. 



4 + 1 



Axi -t- Ax^+i 



< 2- 



A.Ti_i -f Axi 
Axi + Axi+i 



3.3. Revisiting Sweby's diagram. The Sweby's diagram for limiter functions 
[18j are revisited here to include non-uniform grids. The limiter functions considered 
satisfy both the TVD stability condition (3.111 and the order condition (3.21 1. 



For simplicity, it is assumed that the limiter function </'a,_b(') is characterized by 
two parameters A and B which satisfy 



and 



< B < min(2, 2A) 



(3.22) 



(3.23) 



This notation is used extensively in the remaining of the paper to indicate the en- 
hanced limiter functions that take into account of non-uniform grids. With these 
notations the conventional limiters are written as 0(-) = 4'i^i{-), and the correspond- 
ing Sweby's diagram is shown in the Figure |3.5| In this diagram, the admissible 
region of the limiter function is indicated by the shaded region, which is bounded by 
four straight edges: (1) 4>{e) = 2, (2) 0(61) = 29, (3) 0(6') = 1 and (4) 0(6*) = 9. 
The first two edges correspond to the conventional TVD stability conditions given 



by Eqs. (2.13); and the other two edges are obtained from the Lax-Wendroff method 
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,5 2 2.5 3 3.5 




Fig. 3.5. Sweby's diagram for conventional FiG. 3.6. Sweby's diagram for enhanced 

limiter function (f){-) = </>i,i(-) limiter function 4>a,b{') 



[TT] as indicated by Eqs. (3.24) and the Beam- Warming method as indicated by 
Eqs. ( |3.25 ) in the case of uniform meshes in space. 

c< + |c(l-A)«+i-<) 



F 



LW,n 
+ 1/2 



F 



BW.n 
-1/2 



c<+i-ic(l- A) «+!-<) 
c< + ic(l-A)«-<_i) 
c<+i-ic(l-A)«+2-<+i) 



c> 
c < 

c> 
c < 



(3.24) 



(3.25) 



where A = |cAt"|/Ax is the Courant number on the uniform mesh with ceU size Ax. 



In the view of Eq. (3.12) and Eq. (3.8), the two methods given by Eqs. (3.24)-(3.25l 



are equivalent to applying the following "slope limiters" 

1 c> 



0-^^(6') = 



(3.26) 



c < 



c> 



1 c < 



(3.27) 



Thus the admissible region in the conventional Sweby's diagram is explained such that 
the limiter function should be a convex combination of and <b^^ , in addition to 



satisfying Eqs. (2.13) and Eq. (2.16) 



Sweby's diagram is generalized to the limiter function 4)a,b as follows. First, the 
limiter function (f)^^ and 4>^^ are generalized to satisfy Eq. (3.23) as 



B 



c> 
c < 



iBW 



{0) 




c> 
c < 



And the four edges surrounding the admissible region are: (1) 4)a,b{S) — 2, (2) 
0a,b(6') = 261, (3) (t)A,B{0) = B and (4) 0a,s(6') = BO /A; and the corresponding 



Sweby's diagram with the shaded admissible region is illustrated in Figure 3.6 
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Fig. 3.7. Illustration of symmetry-preserving property on uniform meshes: problem 1 (left); 
problem 2 (right) 



Remark: introducing the two parameters A and B that satisfy Eq. (3.22), one 



talks about a class of limiter functions rather than a single one when the name of 
the limiter function is referred to. For example, the enhanced minmod limiter (see 
Section |4]) refers to the class of limiter functions that are generalization to the con- 
ventional minmod limiter function given by Eq. (2.17a) to satisfy Eq. (3.231. In the 
remaining of the paper, the limiter function without subscript like 0™™™°<^ refers 
to the conventional limiter function; whereas the enhanced version is denoted by 



imznmod 



i.e. a class of limiter functions parameterized by A and B. 



3.4. Preserving symmetric solutions. In addition to the TVD stability and 
formal second-order spatial accuracy properties, the slope limiters are also expected 
to preserve the symmetric solutions. On uniform meshes, the condition is given by 



1 



e 



(3.28) 



It is easy to verify that all the limiters listed in Eqs. (2.17) satisfy this condition. 



A geometric interpretation of the symmetry-preserving property of the conven- 
tional limiter function (/)(•) is the following. Consider the REP procedure applied with 
the same uniform mesh to solve two advection problems from time-instance to i^^^ 
for the approximated data in a particular cell fi^ 
1. Problem 1: the governing equation is 



CUr — 



Given data at t^' 



k = 0, ±1, • • • , use limiter function </>(•) to solve for 



the data in Qi at the next time-instance i""*"^, i.e. 
2. Problem 2: the governing equation is 



,n+l 



CVr = 



Given data at t" 



0, ±1, • • • such that 



'i+k 



use the same 



•) to solve for the data in at the next time- instance t"^^, i.e. 



,,n+l 



Then symmetry-preserving is equivalent to u"~^^ — wf^^. This definition is illustrated 



in Figure 3.7 in the case c > 0, in which the nodes represent the cell-averaged values 
at time-instance t" and the slopes represent the reconstructed data after applying the 



limited slopes (c.f. Figure 3.2) 



However, one cannot use the same argument to define the symmetry-preserving 
property when non-uniform grid is concerned, because the mesh is in general non- 
symmetric. To this effect, a definition of symmetry-preserving property for a class of 
slope limiters 0a,s(') parameterized by A and B is proposed as follows 
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II, + = 0. i>I} IV - ml, = 1], t>0 




Fig. 3.8. Illustration of symmetry-preserving property on non-uniform meshes: problem 1 
(left); problem 2 (right) 



Definition 3.2 (Symmetry-preserving). A class of enhanced limiter functions 
0A,b(') 'is said to be symmetry-preserving, if the following two problems give rise to 



^n+l ^ n+1 



1. Problem 1: the governing equation is 

Ut + CUx = 

Given the grid with cell sizes Aa;"^^, and corresponding data at t" : k = 

0, ±1; • • • , solve for u"^^ using the limiter functions (t>A.Bi ) satisfying 

where = 11^+^^^^ <^nd = a.j+Zj^, ' J = * ± 1, ' • ' 

2. Problem 2: the governing equation is 

Vt — CVx = 

Given the grid with cell sizes Aa;^_(_j, and corresponding data at t" : wf^^, , k = 
0, ±1, • • • , such that Aa;^_^j. = and u^j. = """-fe; solve for v^^^ using 

the limiter functions (t>A.B{') satisfying 



where A^ 



. „ , . — and B'" 



J 



I, z ± 1, 



Here the superscripts u or w designate one of the two problems to which a par- 
ticular quantity is relevant. This definition is illustrated in the case c > 0, in analogy 



Eq. (3.8), u 



of F igur e 3.7 for uniform meshes, by Figure 3.8 Let c > 0, using Eqs. (3.7) and 



„n+l 



is equivalent to 



, Ar [(1 - x-mmiu7+i - o - (i - A^i)0^i(ei)K' 

-Ar [(1 - K)<KieT)iv^, - vT) - (1 - A^i)0^i(^^i)(«; 



i+2 



H-l) 
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Employing v'l^^ = <_fc and l^xl^^ = Aa;^.^ for fc = 0, ±1, 



4" 



■nv 



^^1+k 



Aa;'' 



A" 



2— 



i-fe-1 



i — k 



i — k 



^^1+k-l 



"-i+k 



Ax^, + Ax^,+i Aa.^, + Aa;^,_i 

2A2;!' 



Ax" , 



Ax' 



i-fe-l 



^i-k 



It follows that 



1 



Ar [(1 - Ar)c(^r)KVi - o - (i - - ^^i)] 



=<-Ar«-uiLi)+ 



1 



A!' 



(1 - AD^r (^) «-i - <) - (1 - Ar_i)0,iVi (^^) [vju - 



=<-Ar«-u^i)- 



2^^ 



(1 - ^n-^r ( 4) ^rKVi - - (1 - K-i)<t^i+i 



(3.29) 



In the second equality the definition (2.11) is invoked. A sufficient condition for 



Eq. ( |3.29| to hold is the following 

<l^l^k{0) 
9 



Expressing the limiter functions in the form (j)A.B as given in Definition 3.2 this 
condition is equivalent to 



(3.30) 



which suggests the following definition: 

Definition 3.3 (Conjugate limiter function). Given a class of enhanced limiter 
functions (f>A,B characterized by the two parameters A, B as before, the conjugate 
class of limiter functions <j)\ ^ is defined by 



4^*A,B = 4>l/A,B/A 



(3.31) 



is well-defined, since for any A, B that satisfy Eq. (3.22 1, it is always true that 

B min(2, 2 A) 



'<A< A 



min(2, -) 



In addition, it is clear that 
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thus (t>*j^ Q — (j)A.B for any pair A and B. In the view of Eq. (3.30), a sufficient 
condition for the symmetry-preserving property is given by 



'-•AM 



(3.32) 



Finally, in the case of conventional limiter functions (p — the conjugate limiter 
function is 0* = 0| = = 0, thus Eq. (3.32 1 reproduces Eq. (3.28) in the case of 



uniform meshes. 

Remark: an alternative approach to study the symmetry-preserving property is 
given by M. Berger [2 , where a symmetry variable is defined instead of the classical 
smoothness monitors. 



4. Enhanced limiter functions. The conventional limiter function ( 



given 

in Eqs. (2.171 (the superscript name stands for minmod, MC, etc.) is enhanced to a 



class of limiter functions 0™]^'^ so that it falls into the admissible region given by the 
generalized Sweby's diagram (Figure 3.6). The symmetry-preserving condition (3.32) 
is verified after the construction. 



4.1. Enhanced minmod, superbee and MC Hmiters. 



iminmod 



and 



represent lower and upper boundaries of the Sweby's diagram (Figure 3.5). Thus 
direct enhancements are given by 



*-^^^"^o''(0) = max ( 0, J min (6*, A) 



(4.1) 



isuperbee 



(6*) = max ( 0, min {29, B) , min ( 2 



(4.2) 



They are plotted in the generalized Sweby's diagram in Figure 3.6 

The other limiter, 0^*^, is generalized as below (see also Figure (3.6)) 



''Asi^) = max I 0,min I 261 



B 

ATI' 



1),2 



(4.3) 



Supposing 9^0, the symmetry-preserving property of the three generalized 
limiters are easily verified below 



^A,B 



\0) 



iminmod 
^1/A,B/A 



, r. B . ( A 

max I 0, — mm ( 1, — 



^A,B 



, B/A fl 1 
max(0,^mm(^-,^ 



iSuperbee 
'>AB 



— max I 0, mm 12,-1, mm I ^ , ^ 

, 1 BIA\ [B 2 

:max I 0,mm I 2, - • I , mm I — , - 



,'2 B\ fl B/A 
-- max I 0, mm | ^ I , mm I - • ^y^, 2 



I superbee 
~~ ^1/A,BIA 



I superbee, 

'^a'b 
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max 0, min 2 



B 



= max 0, min 



A + 1 
2 B/A (1 



»' \/A + l 



1 ,2 = 



h/A,B/A 



,MC,* 



4.2. Enhanced van Leer limiter. Two generalized van Leer limiters are pro- 
posed in [i2 in terms of symmetry variables, which can be reformulated in the notations 
of this paper as the following 



'A,B 



'{0) 



29 



f) [^r 



A+l 
1 ■ A 



1-1 



< A 



>A 



(4.4) 



iBerger— 2 



(0) 



B(e+i) 

A+l 



Bje+i) 

A+l 



1 - ( 1 - • ^ 

^ ' ^ 8+1 A 



1 _ ( 1 _ 



< A 



> A 



(4.5) 



Both limiters reduce to the conventional van Leer limiter given by Eq. (2.17dl in 



the case A = B = 1; and both satisfy the TVD stability condition, order condition 
and symmetry-preserving condition. They are plotted in Figure [XG} 

A drawback of these two modifications is that they are not smooth at the point 6 = 
A, whereas the conventional 0"'^" ^^^'^ is smooth for all 9 > Q. A smooth enhancement, 
denoted by ^'^'^'^ is provided here. 

Noticing that 



1™ —r, n — 5 — 

k^oc A'' + A'^-i + 



A 



+ A + 1 



= min(l,v4) 



and in the view of Eq. (3.22), there always exists an integer k, such that 



The enhanced van Leer limiter is defined as 



A + l 



(4.6) 



a''+a''-''+--+a+i 



> 



Lvan Leer 



id) 



(4.7) 



gkj^gk-l^ ^e+1 A>'+A''~'^-\ \-A 

9 <Q 

To verify the TVD stability condition, supposing > 0, it is sufhcient to apply 



Eq. (4.6) to show that 



I van Leer ^ 2{9 



qfe-1 



+ 6'fe-i + --- + + l 



< 2min(l,( 



The order condition given by Eq. (3.23) is verified by the following equality 
B{A'' + A'^-i + • ■ 



ivan Leer 
^A,B 



{A) = 



A) A" + A''-^ + --- + A + 1 



+ yl + l A^-" + A^-^ 



A 



= B 
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Finally, the symmetry-preserving condition is verified as following 

^van Leer(g) ^ ^ gk-2 ^ . . . ^ g ^ 1) + 1 + . . . + A + 1 

e ^ 0'' + 0''-^ + --- + 9 + 1 ' Af" + Af"-^ + ■ ■■ + A 
_{B/A)-A- 0''{9-^ + 0-^ + ■■■ + 6*1-'= + 0-'') A''{l + A-^ + --- + A~'') 



6ife(l + + • • • + 0-'') + A-2 + . . . + A-f") 

{B/A){0-^ + 61-2 + ... + e^-^ + 6'-'^) 1 + + • • • + A-*^ 
" l + 6'-i + --- + e'-fe yl-i + A-2 + . . . + yl-/^ 



luari Leer / "'^ 1 J,""" Leer,* I 1 



A sample curve calculated with A; = 3 is plotted in Figure |3.6| And clearly, if A 



B = 1, Eq. (4.6) is satisfied by A: = 1, in which case the enhanced limiter Eq. (4.71 
reduces to the conventional van Leer limiter. 

Remark: This particular limiter can be used as a good example to show that, the 
issue raised in this paper cannot always be resolved by choosing a different smoothness 
monitor and at the same time applying the conventional limiter function. See the proof 
in Appendix [Bj 

4.3. Enhanced van Albada limiter. The van Albada hmiter 0"°" Aibada jg 
sometimes preferred over 0"°" since the former is less restrictive in certain regions 
of 0. The enhanced van Albada limiter proposed here remains a smooth function of 
the smoothness monitor, and is based on the following inequality 

(fc-l)''-! 
fe— s-oo K 



SO that given A and B satisfying Eq. (3.22), there always exists k such that 
Then the limiter function is defined as 



B<2[1+ '—^ ) min(A, 1) (4.8) 



±van Albada f n\ B{6 ' -\- 0^ . . 

To verify the TVD stability condition, one need to show (j)"^^ {0) < 2min(l,6') 
for all 6* > 0, as the following 

The last inequality always holds due to the inequality of geometric and arithmetic 
means 

^-i]0^+ik-i).^^>k(^-iy(^)^0>0 

B J ^ ' (fc-l)S- \B ) \k-lj 



where the last inequality holds because of Eq. (4.8). 
In addition 

B{0^ + 0) ^ f2A \ 1 2 1 f2A \ 1 , 2 1 
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and following the same reasoning one concludes '^"'°''°(6') < 26. 

The order condition 4>^a^ Aibada^^-^ = B is clearly satisfied according to Eq. (4.9) 
and the symmetry-preserving condition is shown as follows 



Ivan Albada 



(0) B{e''-^ + l) {B/A)9''{e-^ +9-'') {B/A){0-'' + 9-^) 



9''+ A 9''{l/A + 9-'') 9-^ + 1/ A 



ivan Albada I \ M""-""- Albada,- 

H/A,B/A \n ] — VA.B 



An example curve with fc = 2 is provided in Figure 3.6 In addition, when fc = 1 



the enhanced limiter Eq. (4.9) coincides with the conventional one given also that 
A = B^l. 



Remark: in practice, evaluating the right hand side of Eq. ( 4.8 ) can be expensive, 
and an economical alternative is due to the observation that 



thus in the numerical examples to be presented, the value k is computed such that 
the following condition is satisfied 



B<2\l + -j min(l,v4) 
5. Numerical examples. 

5.1. Convergence tests. For all the numerical examples here to show the con- 
vergence behavior when the solutions are known to be smooth, the five enhanced 
limiters given in previous section as well as the two modified van Leer limiters by M. 
Berger are tested. Two parameters r = 0.2 and r = 0.3 are used to generate the 
non- uniform meshes as described in Section [21 

5.1.1. One-dimensional smooth Euler problem. The same problem de- 
scribed in Section [2] is solved here with the enhanced limiter functions. With each 
value of the parameter r, five meshes with number of cells ranging from 100 to 1600 
are generated independently to test the convergence behavior. The Li norms of er- 



rors in computed pressure are plotted in Figure 5.1 for the parameters r — 0.2 and 
r = 0.3. Together with the first-order Godunov scheme where no linear reconstruc- 
tion is applied (referred to in the figure as </> = 0), the limiter functions chosen to 
be reflected in the flgures are: conventional minmod limiter, conventional van Albada 
limiter, enhanced minmod limiter, enhanced van Albada limiter and the two modified 
van Leer limiters by M. Berger. 

Further more, all the computed convergence rates are reported in Table |5.1| and 
Table [5^ for the parameter r = 0.2 and r = 0.3 respectively. 



The conclusions from the results are 
1. On non- uniform meshes, the solutions obtained from applying conventional 
limiters, though only show first order convergence, are in general more ac- 
curate than first-order schemes, as shown in Figure |5.1[ MUSCL with the 
enhanced limiters generally show higher-order convergence rates and smaller 
errors than the first order scheme and MUSCL with conventional limiter func- 
tions. 
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Fig. 5.1. Convergence rates in pressure achieved with selected limiter functions on non-uniform 
meshes generated by: r = 0.2 (left) and r = 0.3 (right) 



Table 5.1 

Convergence rates achieved on non-uniform meshes with r = 0.2 by applying enhanced limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


Berger-1 


Berger-2 


p 


1.6603 


1.3193 


2.0490 


2.0935 


1.9620 


2.0888 


2.0742 


u 


1.7309 


1.5222 


2.0645 


2.1103 


2.0645 


2.1063 


2.0978 


p 


1.7308 


1.4729 


2.0566 


2.1119 


2.0023 


2.1088 


2.1003 



Among the enhanced hmiters, (j)™^™'"^ and 0™^'^'"''^ are less accurate than 

the others, thought convergence rates above first-order are observed. 

The enhanced hmiters 4>^%, (/jab^""'', ^vanAibada together with the two 



Umiters (t>^'^g 



iBerger- 
^A.B 



2 



recover second-order convergence in all fluid 



state variables almost perfectly. 



5.1.2. Two-dimensional smooth Euler problem. The two-dimensional test 
is given by the isentropic vortex advection problem |16j . The governing equation is 
the 2D Euler equations 



Wf + f(w)a; + g(w),y = 

f(w)- 



p 

pu 
pv 
E 



pu 
pv? + p 

puv 
u{E+p) 



g(w) 



pv 
puv 
pv^ -\- p 
v{E + p) 



(5.1) 



where u and v are velocity components in x— and y— directions, and the total energy 
per unit volume is hence defined as 



The equation of state is the same as that given in Section [2] for ID Euler equations. 
The computational domain is {x, y) £ [—5, 5] x [—5, 5] with periodic boundary condi- 
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Table 5.2 

Convergence rates achieved on non-uniform meshes with r = 0.3 by applying enhanced limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


Berger-1 


Berger-2 


p 


1.6331 


1.3647 


2.0566 


2.0118 


1.9770 


2.0603 


1.8940 


u 


1.7026 


1.5312 


2.0708 


2.0834 


2.0496 


2.0931 


2.0534 


p 


1.7069 


1.5222 


2.0646 


2.0620 


1.9821 


2.0865 


1.9080 



T-1 



tions. The initial condition is given by 

u{x,y,0) = 1 - ^exp Q(l - - y^) 
v{x,y,0) = 1 + ^exp Q(l - - y^) 

p{x, y, 0)^(l- ^^^^^ exp (1 - - y'] 

pix, y, 0) = (^1 - ^^^^11'^ exp {1-x^- y^] 

and the parameter e is chosen to be e = 5. The initial condition gives rise to an 
isentropic vortex with uniform entropy s = p/p^ = 1; and the vortex's motion is 
described by advection in the (1, 1) direction. After one period of time T — 10.0, the 
exact solution is the same as the initial condition, which is employed to calculate the 
errors of the numerical solutions. 

The problem is solved numerically on non-uniform Cartesian grids, which are 
generated in a way such that the non-uniform spacing in x— and y— directions are 
performed independently according to the algorithm given in Section[2]for the ID case. 
When linear reconstruction is applied, MUSCL technique with chosen slope limiters 
are performed via a dimension-splitting manner. For all the simulations presented, 
second-order TVD Runge-Kutta time-integrator is employed. 

The initial pressure contour is provided in Figure |5.2| together with the solution 
computed using first order Godunov scheme (Figure |5.3[ ); with a non-uniform grid 
generated with the parameter r = 0.2, t he s olutions computed using MUSCL with 
conventional van Albada limiter (Figure 5.4 1 and the enhanced version (p^^-^Aibada 



(Figure 5.5) 



It is evident from the Figures [5.2f|5.5| that, MUSCL technique gives much more 
accurate solutions to the first order scheme; and among the two simulations calculated 
with MUSCL technique, the one using the enhanced van Albada limiter is even more 
accurate than the one using the conventional van Albada limiter. 

Similar as in the previous ID tests, the convergence behavior of Li norms of errors 



in pressure are plotted in Figure 5.6 for selected slope limiters. For all the simulations, 
four meshes with number of cells ranging from 20 x 20 to 160 x 160 are employed for 
the purpose of evaluating the convergence behavior. 

The convergence rates calculated using the two finest meshes are reported in 
Table [573| and Table [5^ for the parameters r = 0.2 and r = 0.3 respectively. 

The following conclusions can be drawn from the numerical results 
1. On non- uniform Cartesian meshes, MUSCL with conventional limiters are 
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-3-2-101234 




Fig. 5.2. Initial pressure for the isentropic Fig. 5.3. Pressure at T = W computed 

vortex advection problem by first order Godunov scheme on non-uniform 

mesh. 





12 3 4 



-4 -3 -2 -1 



2 3 4 



Fig. 5.4. Pressure at T = 10 computed by 
MUSCL with on non-uniform mesh. MUSCL with ,f,j'^gAibada non-uniform mesh. 



Fig. 5.5. Pressure at T = 10 computed 

Lvan / 





Reference mesh size (log scale) 



Reference mesn size (iog scale) 



Fig. 5.6. Convergence rates in pressure achieved with selected limiter functions on non-uniform 
meshes generated by: r = 0.2 (left) and r = 0.3 (right) 



more accurate than first order Godunov scheme, whereas they are less accu- 
rate than MUSCL equipped with the enlianced limiters. 
2. In terms of convergence rates, when tlie meshes are sufficiently refined, the for- 
mal first-order Godunov scheme presents convergence rate below first-order; 
the conventional limiters in MUSCL framework give asymptotic first order 
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Table 5.3 

Convergence rates achieved on non-uniform meshes with r = 0.2 by applying enhanced limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


Berger-1 


Berger-2 


p 


1.8798 


1.5014 


1.8667 


2.0177 


2.1506 


2.0093 


2.0251 


u 


1.7014 


1.6180 


1.9089 


1.9058 


1.8651 


1.9254 


1.9165 


V 


1.6776 


1.6496 


1.9617 


1.9452 


1.9058 


1.9620 


1.9514 


p 


1.8982 


1.5557 


1.9255 


2.1539 


2.2957 


2.1474 


2.1681 



Table 5.4 

Convergence rates achieved on non-uniform meshes with r = 0.3 by applying enhanced limiters 





minmod 


superbee 


MC 


van Leer 


van Albada 


Berger-1 


Berger-2 


p 


1.8339 


1.4765 


1.8499 


1.9765 


2.0095 


1.9949 


2.0254 


u 


1.6970 


1.5680 


1.9060 


1.8959 


1.8178 


1.9246 


1.9054 


V 


1.6599 


1.6058 


1.9610 


1.9224 


1.8052 


1.9523 


1.9279 


p 


1.8549 


1.4981 


1.9130 


2.0986 


2.1158 


2.1303 


2.1654 



convergence; and second order convergences is recovered by the enhanced 
hmiters except for 0™'™°'' and 0^^'^'"^'^, in which case convergence rates 
between 1.5 and 2.0 are observed. 

5.2. Suppression of oscillations near discontinuities. An important pur- 
pose of applying limiters is to enhance the stability property of the numerical algo- 
rithm, especially near discontinuities. The ability of the enhanced limiters to suppress 
spurious oscillations on non-uniform grids is illustrated by ID and 2D benchmark 
problems in this section. For all the test problems, numerical results calculated with 
conventional and enhanced minmod limiter and van Albada limiter are presented; 
applications with other limiter functions lead to similar conclusions. 

5.2.1. One-dimensional shock tube problem. The first problem is the clas- 
sical Sod's shock tube problem 17 , which is governed by the ID Euler equations in 
the domain x G [—2. 2], with initial condition given in primitive variables as 

p{x,0) = 1.0 ( p(x,0) = 0.125 

u(a;,0):=0.0 x<0 I it(a;,0)=0.0 a; > (5.2) 

p(x,0) = 1.0 [ Pix,0) = 0.1 

The problem is solved to T = 0.8. The solution is composed of a left-going rarefaction, 
a right-going contact discontinuity and a right-going shock. The density, velocity and 
pressure profiles at terminal time are plotted in the left columns of Figures |5.7f|5.9[ 
and the enlarged local views are provided in the right columns. All the simulations 
are computed with 200 cells in space and irregularity parameter r = 0.3; the CFL 
number calculated with the reference cell size is chosen to be 0.6. It is clear that the 
three non-linear waves are all resolved well for all the schemes tested. 
It is clear from the figures 

1. The MUSCL schemes (with conventional limiters or enhanced limiters) are 
more accurate than the first order Godunov scheme (that corresponds to the 
case = in the figures) . 

2. No spurious oscillations is visually observable for all the MUSCL schemes, 
either with conventional limiters or with enhanced limiters. This falls in line 
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Fig. 5.7. Density at T = 0.8; global profile (left) and local profile (right) 





1.3 -1.25 -1.2 -1.15 -1.1 -1.05 



-0.95 -0.9 -0.85 



Fig. 5.8. Velocity at T = 0.8; global profile (left) and local profile (right) 





Fig. 5.9. Velocity at T = 0.8; global profile (left) and local profile (right) 



with the previous arguments that the TVD stabihty conditions are the same 
for both uniform and non-uniform meshes. 
3. The enhanced hmiters provide more accurate solutions than the conven- 
tional limiters, which is evident by comparing solution curves computed with 

^rrnnmod g^^^ (j)J^nniod ^ ^j^^g^ computed with and (j)^J^^ ^Ibada ^ 



The second example is the double-shock problem |12j which is governed by ID Euler 
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equations on the domain x G [—3,3] with initial condition given in primitive variables 
p{x,0) = 1.0 ( p{x,0) = 2.0 

u{x,0) = 3.0 x<0 I u(a;,0) = 1.0 a; > (5.3) 
p{x,0) = 1.0 [ p{x,0) = 1.0 

The terminal time of simulation is chosen to be T = 0.8. The solution is composed of a 
slow right-going shock wave, a right-going contact discontinuity and a fast right-going 
shock wave. This is a more challenging problem than the previous one, mainly due to 
the fact that the density peak between the contact discontinuity and the fast shock 
wave is difficult to capture. 150 cells generated with parameter r = 0.3 are employed 
for all the simulations, and CFL number computed with the reference cell size is chosen 
to be 0.6. The density, velocity and pressure profiles are given in Figures [5 . 1 0^5 . 1 2} 
with the left columns being global view and the right columns being local details. 
From the figures one can conclude 

1. The MUSCL schemes (with conventional limiters or enhanced limiters) are 
more accurate than the first order Godunov scheme, especially in resolving 
the density peak near a; = 1.7 as shown in the right part of Figure [5.10| 

2. A few oscillations are observed for pressure between the two shocks, whichever 
limiter function (conventional or enhanced) is applied with MUSCL. This 
violation of the TVD stability can be largely explained by the nonlinearity 
of the Euler equations, and the point here is that both conventional and 
enhanced limiters present similar stability behaviors. 

3. The enhanced limiters lead to more accurate solution in density than the 
conventional counterparts do; whereas the difference in velocity and pressure 
are less distinguishable. 

The last ID example in this section is the Woodward-Colella blast-wave problem 
governed by ID Euler equations on the domain x e [0, 1] with piecewise constant 
initial conditions given in primitive variables 

( p{x,0) = 1.0 
< a; < 0.1 I u{x, 0) = 0.0 0.1 < a; < 0.9 
[ pix,0) = 0.01 

0.9 < a; < 1 (5.4) 



The two boundary conditions are given by reflective walls. The problem is solved up 
to T = 0.038. This problem is characterized by large pressure ratios in the initial 
condition, and a collision of shock waves at around t = 0.027; and the final solution 
presents both large ratios in density and in pressure. Numerical results are computed 
with 200 non-uniform cells generated with r = 0.3, and CFL number calculated with 
reference cell-size fixed to be 0.6. The solution curves are provided in Figures |5.13| - 

Em 

From the figures one observes 
1. The first-order scheme gives totally wrong solution between the two density 
peaks near x — 0.65 and x — 0.8, due to the excessive numerical dissipation; 
MUSCL with conventional limiters lead to roughly correct solution in this 
region, but with small oscillations; and the enhanced limiters compute the 
correct shape, and oscillation-free solutions. In addition, between the two 
enhanced limiters shown, 4>^% Aibada ^gg^^jg more accurate than 0™™"°''. 



p{x,0) 


= 1.0 


u[x, 0) 


= 0.0 


p(a:,0) 


= 1000.0 


p(a;,0) 


^ 1.0 


u{x, 0) 


0.0 


p{x,Q) 


= 100.0 
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Reference 

. - + - ^0 





Fig. 5.10. Density at T = 0.8; global profile (left) and local profile (right) 




Fig. 5.11. Velocity at T = 0.8; global profile (left) and local profile (right) 




Fig. 5.12. Velocity at T = 0.8; global profile (left) and local profile (right) 



2. All MUSCL schemes tested produce more accurate velocities than the first- 
order scheme, and these solutions show similar stability property and are 
hardly distinguishable from each other. 

3. Enhanced limiters lead to more accurate pressure at the peak than the con- 
ventional ones; whereas they are all much more accurate than the first-order 
scheme. 

5.2.2. Two-dimensional Euler problem. Two two-dimensional problems are 
tested. Both problems are governed by the 2D Euler equations. The meshes are 
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obtained in the same way as in the previous isentropic vortex advection example. 
The first problem is a Mach 3 wind tunnel with a step |26| , in which the wind tunnel 
is given by the domain {x, y) € [0, 3] x [0, 1], and the step is placed at a; = 0.6 and with 
height 0.2, as shown in Figures [5 . 1 6H5 . 1 8} Wall boundary conditions are specified at 
upper wall of the wind tunnel, lower wall of the wind tunnel before the step, and the 
left- and upper-boundaries of the step. Incoming flow condition 



p(0,y,i) = 1.4, u{0,y,t)^3.0, v{0,y,t) = 0.0, p{0,y,t) = 1.0 



28 



X. ZENG 




Fig. 5.16. Density at T = 4.0 calculated on: uniform mesh with c/)'"'" Albada (upper-left); non- 
uniform mesh (r = 0.3) with first- order Godunov scheme (upper-right); non-uniform mesh (r = 0.3j 
with <j>'""^ Albada. non-uniform mesh (r = 0.3) with (j,^J^"^ Albada 




Fig. 5.17. Mach number at T = 4.0 calculated on: uniform mesh with cf)^'^" Albada (upper-left); 
non-uniform mesh (r = 0.3) with first-order Godunov scheme (upper- right); non-uniform mesh 
(r = 0.3) with cji""" Albada. non-uniform mesh (r = 0.3) with ff,^^n^Albada 



is specified at the left opening of the wind tunnel; whereas the outgoing flow condition 
at a; = 3.0 is given by homogeneous Neumann boundary condition 

p,(3,y,i) = 0.0, M,(3,?/,t) =0.0, w,(3,?/,0 =0.0, p,(3, = 0.0 

At the terminal time T = 4.0, the solution presents shock-shock interactions, shock- 
boundary layer interactions and shock-wall reflections, where the boundary layer is 
generated at the corner of the step which is a singular point in the flow. Numerical 
solutions computed on both uniform and non-uniform (with r = 0.3) 150 x 50 (which 
number includes those immersed in the step) meshes are presented as 50 contours 



for density, Mach number and pressure in Figure 5.16 Figure 5.17 and Figure 5.18 
respectively. Among the four simulations, MUSCL with Aibada applied for 
the uniform mesh, and the results are reported in the upper-left parts of the figures. 
For the non-uniform meshes, results computed from first order Godunov scheme are 
presented in the upper-right parts of the figures; MUSCL with conventional limiter 
^van Albada ^.j^g lower-lcft parts, and MUSCL with enhanced limiter (/)'"/'% are 
in the lower-right parts. 

From the figures one conclude 

1. All the simulations lead to stable results. 

2. The first order scheme produce wrong results, since the interaction between 
the first shock wave and the reflected shock wave from upper wall is not 
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captured at all. 

3. Using the numerical results computed on the uniform mesh (upper-left parts) 
as reference, MUSCL with enhanced limiter Aibada igg^^g almost iden- 
tical contours; whereas MUSCL with conventional limiter 0'""" ^'^f"^" gives 
rise to less accurate solutions, especially in the following regions: 1. After the 
shock front on the lower wall in front of the step (near the location (0.3, 0.0) 
on the lower wall); 2. After the first reflected shock from upper wall (near 
the location (0.6, 1.0) on the upper wall). 
The second problem is the double Mach reflection of a Mach 10 shock The 
physical setup for the problem can be found in the paper cited; the numerical problem 
setup is given below. The computational domain is (x, y) e [0,4] x [0,1], which is 
discretized by uniform or non-uniform Cartesian meshes. The initial condition is given 
by flow separated by a Mach 10 shock, located along the straight line x — l/Q + y/\/?>. 
The flow state vector to the left of the shock is given by the post-shock status 



p(a:,2/,0) = 8.0 
u{x,y,Q) 

?;(x,2;,0) = -8.25 sin (f) 
p{x, y, 0) = 116.5 



.25 cos (l; 



X < 



V 

V3 



(5.5) 



and the flow state vector to the right of the shock is given by the pre-shock status 



pix,y,0) = 1.4 
u{x,y,0) = 0.0 
i;(x,2/,0) =0.0 
p(x,y,0) = 1.0 



X > 



y 

^/3 



(5.6) 



The in-flow condition at a; = and flow condition on the lower boundary before 
a; = 1/6 are given by the post-shock fluid state vector Eqs. (5.5) for all time; the out- 



flow boundary condition at a; = 4 is the homogeneous Neumann boundary condition 

p,(4,y,i) = 0.0, u^(4,y,t)=0.0, v^{A,y,t)^0.0, p,(4, y, i) - 0.0 

On the lower boundary after the point x — 1/6, reflective wall boundary condition 
is imposed; and on the uppe r boundary the exact flow solution is imposed, i.e. the 
post -sho ck values Eqs. (5.5) for x < 1/6 + \/3(l -|- 20t), and the pre-shock values 
Eqs. [Hj] for a; > 1/6 + \/3(l + 20t) . 
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Fig. 5.19. Density at T = 0.2 calculated on: uniform mesh with c/)'"'" Albada (upper-left); non- 
uniform mesh (r = 0.3j with first-order Godunov scheme (upper-right); non-uniform mesh (r = 0.3j 
with 0""" Albada. non-uniform mesh (r = O.Sj with (j,^J^"^ Albada 




Fig. 5.20. Mach number at T = 0.2 calculated on: uniform mesh with cf)^'^'" Albada (uppgj-left) ; 
non-uniform mesh (r = 0.3) with first-order Godunov scheme (upper- right); non-uniform mesh 
(r = 0.3 j with cji""" Albada. non-uniform mesh (r = O.SJ with ^'^ri^ Albada 



The problem is solved to T — 0.2, and the numerical solutions in the domain 



[0,3] X [0,1] obtained from four simulations are presented in Figures 5.19 5.21 For 



all the simulations, the size of the mesh is 480 x 120. The upper-left parts represent the 
solutions obtained from MUSCL with (/i'"™ -^'fcada uniform grid. The other three 
simulations use non-uniform grids generated with r = 0.3: the upper-right parts are 
obtained from first-order Godunov scheme; the lower-left figures are computed using 
MUSCL with the conventional limiter 0'"'" ^""^''<^; the lower-right contours are results 
using MUSCL with the enhanced limiter t/)™^ Albada ^ should be mentioned that the 
spatial resolution used here is not fine enough to resolve the fine structures, especially 
near the wave interaction sections; however the results showing here are sufficient for 
the comparison purpose. 

It can be observed from the figures that 

1. All simulations lead to stable results. 

2. For this more challenging problem, the accuracy is mainly determined by 
the resolution of the meshes rather than the former order of the numerical 
schemes. However, one can still observe that the first scheme is less accurate, 
especially near the propagating shock fronts. 

3. The structure in the shock reaction region computed from the uniform mesh 
is different from those obtained from the non-uniform mesh. This is mainly 
due to the fact that the three simulations use the same non-uniform meshes, 
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which is consistent with the fact that the mesh is not fine enough to resolve 
the physical structure of the flow solutions. 

5.3. Conclusions from the numerical examples. To summarize the numer- 
ical results 

1. When the solution is smooth, enhanced limiters with MUSCL recover second- 
order convergence rates in space whereas conventional limiters only lead to 

first-order space accuracy. However, MUSCL with either enhanced or conven- 
tional limiters leads to more accurate solutions than the first order Godunov 
scheme. 

2. For solutions that exhibit discontinuities, both conventional and enhanced 
limiters show similar stability property, which is consistent with the argu- 
ments before. 

3. For very challenging problems like the double-Mach reflection problem, when 
the fine fiow structure is under-resolved, the inaccuracy mainly comes from 
the mesh resolution rather than the formal order space accuracy of the nu- 
merical scheme. And for such problems, applying conventional limiters or 
enhanced limiters lead to visually similar solutions. 

Especially, the last observation partially explains why it is hard to discover in 
practice that conventional limiters are inappropriate when applied on non-uniform 
meshes. However, this paper shows that by applying the enhanced limiters, which 
only needs minor modifications to the conventional ones, the accuracy in smooth 
regions can be greatly improved. 

6. Conclusions. A general approach of studying and enhancing the limiters 

on non-uniform grids is presented. The theory presented follows the REP approach 
and Harten's stability theory, with modifications taking into account the non-uniform 
computational grids. The analysis is conducted for ID problems, and it explains why 
formal second-order spatial accuracy is lost when conventional limiters are applied 
with non-uniform grids, and at the same time the TVD stability is retained. This 
study provides sufficient conditions for slope limiters on non-uniform meshes, so that 
TVD stability, formal second-order spatial accuracy and symmetry-preserving prop- 
erty are achieved in the resulting numerical schemes. The results extend naturally 
to multiple dimensions, provided that cell-centered finite volume methods arc used. 
Several most widely used limiter functions, including the van Leer and van Albada 
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limiters, are enhanced to satisfy these conditions; in addition, the latter two retain 
their smoothness property in the enhanced form. The enhanced hmiters are apphed to 
solve both ID and 2D benchmark problems to show that second-order spatial accuracy 
is indeed recovered on non-uniform meshes, and TVD stability is retained. 
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Appendix A. Measuring errors and computing convergence rates on 
non-uniform meshes. Given two randomly generated irregular meshes of the do- 
main [a, 6], which are defined by the cell faces 

JVIesh I: a — x{^2 < < ■ ■ • < ^{Ni-i)+i/2 < ^Wi+i/2 — ^ 
Mesh II: a = x{j2 < < ■ ■ < 2;(j(/2_i)+i/2 < 2:^2+1/2 ^ ^ 

The reference cell-sizes for the two meshes are h' = {b- a) /Ni and h" = (5 - a) /N2 
respectively. 

For generic variable u to be solved for, define uf and uj^ as the numerical solutions 
computed by the same scheme on the two meshes respectively at a given terminal 
time T. In addition, define u^'^^{x,T) as the reference solution at location x and 
terminal time T is available either from analytical solution or from numerical solution 
computed on a much finer grids. The Li norm of the computed errors are hence 
calculated as 



1=1 



Then the convergence rates is estimated by 

log(i^^^-bg(i?^ 

\og{hi) - log(/i") ^ ■ ^ 

Note that the equation above is widely used for estimating convergence rates on 
uniform meshes, and its validity for the same purpose on non-uniform meshes is open 
for discussion. However the numerical results presented in Section [5] justifies this 
choice for irregular meshes generated in the way described in Section [2j 

Appendix B. Effect of alternative smoothness monitors on the van Leer 
limiter. An alternative strategy to deal with the issues raised in this paper is to 



choose different smoothness monitors from that given by Eq. (2.11), while still us- 
ing conventional limiter functions. These smoothness monitors could be designed to 
account for non-uniform meshes, like 



Mj— ttj-l A A A A 

^'+^~"' Axi-i + Axi — Ui Axi_i -t- Axi 



(B.l) 



and accordingly the limited slope u" is 



^ " ' x.,+i -X, ' ' Ax, + Ax,+i ^ ' 



34 



X. ZENG 



Using these alternatives with the conventional minmod, superbee and MC limiters, 
the effect is the same as the enhancement proposed in this paper. On the other hand, 
the same does not necessarily hold for general slope limiters. For example, see the 
following result concerning the van Leer limiter and a large range of smoothness 
monitors. 

Theorem B.l. Supposing any alternative smoothness monitor and limited slopes 
given by the following 

Of = a^ef (B.3) 

(B.4) 



ivan Leer ( nalt 



Aa 



where 6i is defined by Eq. \2.11 ). a; and (3i are positive constants that only rely on 
the grid information; and p Cz {±1}, q€ {0,1}. Eqs. (B.3 B.4) are equivalent to 



I van Leer i+1 



Aa 



whe 



ivan Leer 



van heer 



(B.5) 



Then the equivalent limiter 0™" ^'^'^'^ cannot satisfy both the order condition (3.21) 



and the TVD stability condition [3.11) for arbitrary non-uniform grids. 

Before going to the proof of this theorem, the two parameters p and q is motivated 
from the following 



B.3 


B.4 


1, Theorem 


B.l 



possible smoothness monitors and reconstructed slopes 

Proof. [Proof of Theorem 
equivalent limiter given by Eq. 



B.l 



Suppose 6i > in the following proof. The 



B.5| has the following explicit expression 



ivan Leer 



1 + Qfi 



Noting p ^ ±1 and q = 0,1, the formula above can be simplified as 

nt 

Avan Leer ( -i o 

a + 06' 

where a and h are positive constants that are determined by a^, Pi and p. Define the 



grid parameters A 



Thus it is sufficient to proof 



^"'-^+^"' and B = . ^^f^' 
that one cannot choose a and b as functions of A and B, such that the following two 
conditions hold for all 9 > Q, and all possible pair {A,B). 



< 



< 2min(l, 



(B.6) 



and also 



A* 



bA 



B 



(B.7) 
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First consider Eq. (B.6I, one obtains 
at 



a + b0 



< 2, > ^ t = 1, 26 > 1 



Applying t = 1, it follows again from Eq. (B.6I that 



a + be 



< 29, V6' > 



2a > 1 



Next consider Eq. (B.7), by applying i = 1, it leads to 

2A 



2a + 2bA 



B 



Thus if both Eq. (B.6I and Eq. (B.7 1 hold, one must have 



2A 
I3 



>l + A, VA,S, s.t. < B < 2min(l,y4) 



which is, however, not true. □ 



